home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / brent_cc.lha / brent_cc / fminbr.cc < prev    next >
C/C++ Source or Header  |  1993-08-08  |  6KB  |  166 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *              Numerical Math Package
  6.  *
  7.  *            Brent's one-dimensional minimizer 
  8.  *
  9.  *         finds a local minimum of a single argument function
  10.  *              over the given range
  11.  *
  12.  * Input
  13.  *    double fminbr(ax,bx,f,tol)
  14.  *    const double ax                a and b, a < b, specify the interval
  15.  *    const double bx          the minimum is to be sought in
  16.  *    double (*f)(const double x)    Ptr to the function under investigation
  17.  *    const double tol        Acceptable tolerance for the minimum
  18.  *                    location. It is an optional parameter
  19.  *                    with default value EPSILON
  20.  *
  21.  * Output
  22.  *    Fminbr returns an estimate to the minimum location with accuracy
  23.  *    3*SQRT_EPSILON*abs(x) + tol.
  24.  *    The procedure always determines a local minimum, which coincides with
  25.  *    the global one if and only if the function under investigation is
  26.  *    unimodular.
  27.  *    If a function being examined possesses no local minimum within
  28.  *    the given range, Fminbr returns either the left or the right end
  29.  *    point of the interval, wherever the function value is smaller.
  30.  *
  31.  * Algorithm
  32.  *    G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
  33.  *    computations. M., Mir, 1980, p.202 of the Russian edition
  34.  *
  35.  * The function makes use of the "gold section" procedure combined with
  36.  * the parabolic interpolation.
  37.  * At every step program operates three abscissae - x,v, and w.
  38.  *     x - the last and the best approximation to the minimum location,
  39.  *        i.e. f(x) <= f(a) or/and f(x) <= f(b)
  40.  *         (if the function f has a local minimum in (a,b), then both
  41.  *           conditions are met after one or two steps).
  42.  *    v,w are previous approximations to the minimum location. They may
  43.  *    coincide with a, b, or x (although the algorithm tries to make all
  44.  *    u, v, and w distinct). 
  45.  * Points x, v, and w are used to construct an interpolating parabola,
  46.  * whose minimum is to be regarded as a new approximation to the minimum
  47.  * location if the former falls within [a,b] and reduces the minimum 
  48.  * encompassing interval to a larger extent than the gold section procedure.
  49.  * When f(x) has a second derivative positive at the point of minimum
  50.  * (not coinciding with a or b) the procedure converges superlinearly
  51.  * at a rate of about 1.324
  52.  *
  53.  ************************************************************************
  54.  */
  55.  
  56. #include "math_num.h"
  57.  
  58. #pragma implementation "math_num.h"
  59.  
  60.  
  61. double fminbr(                // An estimate to the min location
  62.     const double ax,        // Specify the interval the minimum
  63.     const double bx,        // to be sought in
  64.     double (*f)(const double x),    // Function under investigation
  65.     const double tol = EPSILON )    // Acceptable tolerance
  66. {
  67.   double a = ax, b = bx;        // Current interval
  68.   double x,v,w;                // Abscissae, descr. see above
  69.   double fx;                // f(x)
  70.   double fv;                // f(v)
  71.   double fw;                // f(w)
  72.   const double r = (3-sqrt(5.0))/2;    // Gold section ratio
  73.  
  74.   assure( tol > 0, "Tolerance must be positive");
  75.   assure( b > a, 
  76.      "Left end point of the interval should be strictly less than the "
  77.      "right one" );
  78.  
  79.   v = a + r*(b-a);  fv = (*f)(v);       // First step - always gold section
  80.   x = v;  w = v;
  81.   fx=fv;  fw=fv;
  82.  
  83.   for(;;)        // Main iteration loop
  84.   {
  85.     const double range = b-a;        // Interval over which the minimum
  86.                     // is sought for
  87.     const double midpoint = (a+b)/2;
  88.     const double tol_act =        // Actual tolerance
  89.         SQRT_EPSILON*fabs(x) + tol/3;
  90.  
  91.        
  92.  
  93.     if( fabs(x-midpoint) + range/2 <= 2*tol_act )
  94.       return x;                // Acceptable approximation is found
  95.  
  96.                     // Compute a new step with the gold
  97.                     // section
  98.     double new_step = r * ( x < midpoint ? b-x : a-x );
  99.  
  100.  
  101.                 // Decide on the interpolation  
  102.     if( fabs(x-w) >= tol_act  )        // If x and w are distinct
  103.     {                    // interpolatiom may be tried
  104.       register double p;         // Interpolation step is calcula-
  105.       register double q;                  // ted as p/q; division operation
  106.                                         // is delayed until last moment
  107.       register double t;
  108.  
  109.       t = (x-w) * (fx-fv);
  110.       q = (x-v) * (fx-fw);
  111.       p = (x-v)*q - (x-w)*t;
  112.       q = 2*(q-t);
  113.  
  114.       if( q > 0 )            // Formulas above computed new_step
  115.     p = -p;                // = p/q with wrong sign (on purpose).
  116.       else                // Correct this, but in such a way so
  117.     q = -q;                // that q would be positive
  118.  
  119.       if( fabs(p) < fabs(new_step*q) &&    // If x+p/q falls in [a,b] and is not
  120.      p > q*(a-x+2*tol_act) &&    // too close to a and b, and isn't
  121.      p < q*(b-x-2*tol_act)  )       // too large, it is accepted
  122.        new_step = p/q;
  123.                     // If p/q is too large then the
  124.                     // gold section procedure would
  125.                     // reduce [a,b] to larger extent
  126.     }
  127.  
  128.     if( fabs(new_step) < tol_act )    // Adjust the step to be not less
  129.       if( new_step > 0 )        // than tolerance
  130.     new_step = tol_act;
  131.       else
  132.     new_step = -tol_act;
  133.  
  134.                 // Obtain the next approximation to min
  135.                     // and reduce the encompassing interval
  136.     register double t = x + new_step;    // Tentative point for the min
  137.     register double ft = (*f)(t);
  138.     if( ft <= fx )
  139.     {                                     // t is a better approximation
  140.       if( t < x )            // Reduce the interval so that
  141.     b = x;                            // t would fall within it
  142.       else
  143.     a = x;
  144.       
  145.       v = w;  w = x;  x = t;        // Assign the best approx to x
  146.       fv=fw;  fw=fx;  fx=ft;
  147.     }
  148.     else                                  // x remains the better approx
  149.     {
  150.       if( t < x )            // Reduce the interval encompassing x
  151.     a = t;                   
  152.       else
  153.     b = t;
  154.       
  155.       if( ft <= fw || w==x )
  156.       {
  157.     v = w;  w = t;
  158.     fv=fw;  fw=ft;
  159.       }
  160.       else if( ft<=fv || v==x || v==w )
  161.     v = t, fv = ft;
  162.     }
  163.   }        // ===== End of loop =====
  164.  
  165. }
  166.